#include "header.h"

/* Corrects density and pressure after NWBC  */

/* ------------------------------------------------------------------- */
void correct_nwbc(int point,struct node* p,double** var_old){
/* ------------------------------------------------------------------- */
int i,j,neigh;
double dx,dy,ds2,weight,var[4];

double H,S,So,po,Ho,u,v,w;
double num,denom;
double Ssource,psource,Hsource,cost,sint;
double Sold,Pold,Hold,Srelax,Prelax,Hrelax;
double sx, sy, sz, sxp, syp, szp;

Srelax = Prelax = Hrelax = swbc_relax;

Pold  = var_old[point][4];

sx = sy = sz = 0.0;
sxp = syp = szp = 0.0;

for(i=0;i<p[point].total_conn;i++){
	neigh = p[point].Conn[i];

	sx  += p[point].X[i];
	sy  += p[point].Y[i];
	sz  += p[point].Z[i];

	sxp += p[point].X[i]*p[neigh].var[4];
	syp += p[point].Y[i]*p[neigh].var[4];
	szp += p[point].Z[i]*p[neigh].var[4];
    }

/* --------------------- Strong boundary condition -------------- */
u    = var_old[point][1];
v    = var_old[point][2];
w    = var_old[point][3];
denom = p[point].n[0]*sx + p[point].n[1]*sy + p[point].n[2]*sz;
if(denom == 0.0){
    printf("solid: Denominator zero at node %d\n",point);
    exit(0);
    }
/*--------------------- Normal Momentum equation ---------------------*/
psource = var_old[point][0]*( p[point].n[0]*( u*p[point].qx[1] + v*p[point].qy[1] + w*p[point].qz[1])
                            + p[point].n[1]*( u*p[point].qx[2] + v*p[point].qy[2] + w*p[point].qz[2])
                            + p[point].n[2]*( u*p[point].qx[3] + v*p[point].qy[3] + w*p[point].qz[3]));
/*
psource = var_old[point][0]*( cost*( u*p[point].qx[1] + v*p[point].qy[1])*p[point].phi[1]
                            + sint*( u*p[point].qx[2] + v*p[point].qy[2])*p[point].phi[2]);
							*/
num = p[point].n[0]*sxp + p[point].n[1]*syp + p[point].n[2]*szp + psource;
po  = num/denom;
po  = Pold + Prelax*(po - Pold);
if(po <= 0.0){
    printf("correct_nwbc: Pressure is negative at node %d\n",point);
    printf("pressure = %f\n",po);
    printxy(point);
    exit(0);
    }
/*--------------------- Enthalpy equation ---------------------*/
p[point].var[4] = po;

}
